# This script will demonstrate how to analyze the amphibian data collected as 
# part of the LongTerm Ecological Monitoring Initiative

#
# Only one study area at time can only be analyzed with a script. 
#
# This was programmed by Carl James Schwarz, Statistics and Actuarial Science, SFU
# cschwarz@stat.sfu.ca
#
# 2017-03-27 Separated egg mass by species. Imputed 0 counts when species not seen in a year
# 2017-02-28 First Edition

# Summary of calling protocol
#    Define a survey transect along or through your wetland. … 
#    Surveyors visit the monitoring site(s) in spring and listen for calling males, 
#    recording the species and approximate number of each. 
#    Repeat surveys increase the probability that species will be detected.”

# Summary of visual protocol
#    Surveyors monitor breeding site(s) during the active season (spring to fall),
#    walking the shoreline of a wetland recording all species and life stages 
#    encountered. These include egg masses, tadpoles and adults. 
#    Repeat surveys increase the probability that species will be detected.”

# A separate analysis is done on EACH species present


# load libraries
library(car)       # for testing for autocorrelation (2 libraries needed - see dwtest)
library(ggfortify) # for residual and other diagnostic plot
## Loading required package: ggplot2
library(ggplot2)   # for plotting
library(lmtest)    # for testing for autocorrelation
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(plyr)      # for group processing
library(readxl)    # for opening the Excel spreadsheets and reading off them
library(reshape2)  # for melting and casting
library(lmerTest)  # for the linear mixed modelling
## Loading required package: Matrix
## Loading required package: lme4
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step
library(stringr)   # string handling (like case conversion)

# Load some common functions
source("../CommonFiles/common.functions.R")


cat("\n\n ***** Amphibian Analysis - Single Site Analysis *****  \n\n")
## 
## 
##  ***** Amphibian Analysis - Single Site Analysis *****
# get the data from the Excel work.books.
# we put the list of work books here, including the file type (xls or xlsx).
# You can put multiple stations here because the station information is included on the raw data

work.books.csv <- textConnection(
"file.name
General Survey Using Transects-amphibians_AliceLake2013.xls
General Survey Using Transects-amphibians_AliceLake2014.xls
General Survey Using Transects-amphibians_AliceLake2015.xls
")

work.books <- read.csv(work.books.csv, as.is=TRUE, strip.white=TRUE, header=TRUE)
cat("File names with the data \n")
## File names with the data
work.books
##                                                     file.name
## 1 General Survey Using Transects-amphibians_AliceLake2013.xls
## 2 General Survey Using Transects-amphibians_AliceLake2014.xls
## 3 General Survey Using Transects-amphibians_AliceLake2015.xls
# read each workbook and put all of the data together into one big data frame
amph.df <- plyr::ddply(work.books, "file.name", function(x){
  file.name <- file.path("Data",x$file.name)
  data <- readxl::read_excel(file.name, sheet="General Survey")
  data
})
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 0d 3b 00 00 00 00 f8 00 00 00 04 00 
## DEFINEDNAME: 21 00 00 01 0b 00 00 00 06 00 00 00 00 00 00 0d 3b 01 00 00 00 00 00 00 00 04 00 
## DEFINEDNAME: 20 00 00 01 0b 00 00 00 05 00 00 00 00 00 00 07 3b 00 00 00 00 00 00 00 00 ff 00
# fix up variable names in the data.frame.
# Variable names in R must start with a letter and contain letters or number or _. 
# Blanks in variable names are not normally allowed. Blanks will be replaced by . (period)
cat("\nOriginal variable names in data frame\n")
## 
## Original variable names in data frame
names(amph.df)
##  [1] "file.name"                              
##  [2] "Study Area Name"                        
##  [3] "Transect Label"                         
##  [4] "Date"                                   
##  [5] "Time"                                   
##  [6] "End Time"                               
##  [7] "Surveyor"                               
##  [8] "Species"                                
##  [9] "Count"                                  
## [10] "UTM Zone"                               
## [11] "Easting"                                
## [12] "Northing"                               
## [13] "Survey Type"                            
## [14] "Life Stage"                             
## [15] "Background Noise (for auditory surveys)"
## [16] "Calling Noise (for auditory surveys)"   
## [17] "Preceeding 24 hour precipitation"       
## [18] "Air temp (c)"                           
## [19] "Cloud cover"                            
## [20] "Wind speed"                             
## [21] "Snow coverage on ground"                
## [22] "Current precipitation"                  
## [23] "BEC"                                    
## [24] "Comments"                               
## [25] "Preceeding 24 hour precipitation (mm)"
names(amph.df) <- make.names(names(amph.df))

cat("\nCorrected variable names of data frame\n")
## 
## Corrected variable names of data frame
names(amph.df)
##  [1] "file.name"                              
##  [2] "Study.Area.Name"                        
##  [3] "Transect.Label"                         
##  [4] "Date"                                   
##  [5] "Time"                                   
##  [6] "End.Time"                               
##  [7] "Surveyor"                               
##  [8] "Species"                                
##  [9] "Count"                                  
## [10] "UTM.Zone"                               
## [11] "Easting"                                
## [12] "Northing"                               
## [13] "Survey.Type"                            
## [14] "Life.Stage"                             
## [15] "Background.Noise..for.auditory.surveys."
## [16] "Calling.Noise..for.auditory.surveys."   
## [17] "Preceeding.24.hour.precipitation"       
## [18] "Air.temp..c."                           
## [19] "Cloud.cover"                            
## [20] "Wind.speed"                             
## [21] "Snow.coverage.on.ground"                
## [22] "Current.precipitation"                  
## [23] "BEC"                                    
## [24] "Comments"                               
## [25] "Preceeding.24.hour.precipitation..mm."
# Convert dates to R date format
xtabs(~Date, data=amph.df, exclude=NULL, na.action=na.pass)  # check the date formats. Make sure that all yyyy-mm-dd
## Date
## 2013-05-22 2013-06-10 2014-04-09 2015-04-10 
##         11          6         12          9
amph.df$Date <- as.Date(amph.df$Date, "%d-%b-%y", tz="UTC")
amph.df$Year <- as.numeric(format(amph.df$Date, "%Y"))



# Check that the Study Area Name is the same across all years
# Look at the output from the xtabs() to see if there are multiple spellings 
# of the same Study.Area.Name.

# We will convert the Study.Area.Name to Proper Case.
amph.df$Study.Area.Name <- stringr::str_to_title(amph.df$Study.Area.Name)
xtabs(~Study.Area.Name,      data=amph.df, exclude=NULL, na.action=na.pass)
## Study.Area.Name
## Alice Lake Park 
##              38
xtabs(~Study.Area.Name+Year, data=amph.df, exclude=NULL, na.action=na.pass)
##                  Year
## Study.Area.Name   2013 2014 2015
##   Alice Lake Park   17   12    9
# Check the Transect.Labels for typos
xtabs(~Study.Area.Name+Transect.Label, data=amph.df, exclude=NULL, na.action=na.pass)
##                  Transect.Label
## Study.Area.Name   Fawn Lake Station Fawn Lake Transect
##   Alice Lake Park                 3                 35
xtabs(~Transect.Label+Year,            data=amph.df, exclude=NULL, na.action=na.pass)
##                     Year
## Transect.Label       2013 2014 2015
##   Fawn Lake Station     2    1    0
##   Fawn Lake Transect   15   11    9
# When is each station measured?
xtabs(~Transect.Label+Year, data=amph.df, exclude=NULL, na.action=na.pass)
##                     Year
## Transect.Label       2013 2014 2015
##   Fawn Lake Station     2    1    0
##   Fawn Lake Transect   15   11    9
xtabs(~Date+Transect.Label, data=amph.df, exclude=NULL, na.action=na.pass)
##             Transect.Label
## Date         Fawn Lake Station Fawn Lake Transect
##   2013-05-22                 1                 10
##   2013-06-10                 1                  5
##   2014-04-09                 1                 11
##   2015-04-10                 0                  9
# Check the Species code to make sure that they are all ok
# This isn't used anywhere in the analysis but is useful to know
xtabs(~Species, data=amph.df, exclude=NULL, na.action=na.pass)
## Species
##      A-AMGR      A-AMMA      A-ANBO      A-PRSE      A-PSRE      A-RAAU 
##           5           2           6           8           9           5 
##   AMPHIBION Unidentifed 
##           2           1
xtabs(~Species+Year, data=amph.df, exclude=NULL, na.action=na.pass)
##              Year
## Species       2013 2014 2015
##   A-AMGR         4    0    1
##   A-AMMA         2    0    0
##   A-ANBO         1    1    4
##   A-PRSE         0    8    0
##   A-PSRE         6    1    2
##   A-RAAU         2    2    1
##   AMPHIBION      2    0    0
##   Unidentifed    0    0    1
# Check the survey type fiels
xtabs(~Year+Survey.Type,  data=amph.df, exclude=NULL, na.action=na.pass)
##       Survey.Type
## Year   AU VI
##   2013  2 15
##   2014  1 11
##   2015  0  9
# Check the Life Stage Code
xtabs(~Year+Life.Stage,  data=amph.df, exclude=NULL, na.action=na.pass)
##       Life.Stage
## Year   AD EG LA
##   2013  3  8  6
##   2014  3  9  0
##   2015  4  4  1
xtabs(~Survey.Type+Life.Stage, data=amph.df, exclude=NULL, na.action=na.pass)
##            Life.Stage
## Survey.Type AD EG LA
##          AU  3  0  0
##          VI  7 21  7
# list the adult records
amph.df[ amph.df$Life.Stage=="AD",c("Study.Area.Name","Date","Survey.Type","Life.Stage","Count")]
##    Study.Area.Name       Date Survey.Type Life.Stage Count
## 1  Alice Lake Park 2013-05-22          AU         AD     3
## 7  Alice Lake Park 2013-05-22          VI         AD     2
## 17 Alice Lake Park 2013-06-10          AU         AD     1
## 18 Alice Lake Park 2014-04-09          AU         AD    NA
## 19 Alice Lake Park 2014-04-09          VI         AD     3
## 26 Alice Lake Park 2014-04-09          VI         AD     1
## 33 Alice Lake Park 2015-04-10          VI         AD    10
## 34 Alice Lake Park 2015-04-10          VI         AD     5
## 36 Alice Lake Park 2015-04-10          VI         AD     1
## 38 Alice Lake Park 2015-04-10          VI         AD     1
# Get the file prefix
file.prefix <- make.names(amph.df$Study.Area.Name[1])
file.prefix <- gsub(".", '-', file.prefix, fixed=TRUE) # convert blanks to -
file.prefix <- file.path("Plots", file.prefix)
#  Analysis of the calliing data

# Not possible at this time. See document.
#  Analysis of the visual survey protocol

# Look at # egg masses 
# We need to assume that detection probabilities are consistent over time and space
# There is only one monitoring station at Alice Lake so we don't need to model transect or year-specific
# process error. But I've generalized the code so that it also works with more than one station.

# Must impute 0 values for species not present in a year

dim(amph.df)
## [1] 38 26
amph.v.df <- amph.df[ amph.df$Survey.Type=='VI',]
dim(amph.v.df)
## [1] 35 26
# Count the total number of egg masses seen by species
eggmass.count.by.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
                       n.eggmass=sum(Life.Stage=='EG'))
eggmass.count.by.species
##    Study.Area.Name     Transect.Label Year     Species n.eggmass
## 1  Alice Lake Park Fawn Lake Transect 2013      A-AMGR         3
## 2  Alice Lake Park Fawn Lake Transect 2013      A-AMMA         2
## 3  Alice Lake Park Fawn Lake Transect 2013      A-ANBO         0
## 4  Alice Lake Park Fawn Lake Transect 2013      A-PSRE         2
## 5  Alice Lake Park Fawn Lake Transect 2013      A-RAAU         1
## 6  Alice Lake Park Fawn Lake Transect 2013   AMPHIBION         0
## 7  Alice Lake Park Fawn Lake Transect 2014      A-ANBO         0
## 8  Alice Lake Park Fawn Lake Transect 2014      A-PRSE         7
## 9  Alice Lake Park Fawn Lake Transect 2014      A-RAAU         2
## 10 Alice Lake Park Fawn Lake Transect 2015      A-AMGR         1
## 11 Alice Lake Park Fawn Lake Transect 2015      A-ANBO         1
## 12 Alice Lake Park Fawn Lake Transect 2015      A-PSRE         1
## 13 Alice Lake Park Fawn Lake Transect 2015      A-RAAU         1
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed         0
# Count the total number of egg masses seen for ALL species
eggmass.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
                       n.eggmass=sum(Life.Stage=='EG'))
eggmass.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year"), summarize,
                       n.eggmass=sum(Life.Stage=='EG'))
eggmass.count.all.species$Species <- 'ALL.species'

eggmass.count <- rbind(eggmass.count.by.species, eggmass.count.all.species)




# Need to impute 0 values for species not seen in a year


# check out records where eggmass count not identified
select <- eggmass.count$Species %in% c("AMPHIBION","Unidentifed") # Notice type in unidentified
eggmass.count[ select,]
##    Study.Area.Name     Transect.Label Year     Species n.eggmass
## 6  Alice Lake Park Fawn Lake Transect 2013   AMPHIBION         0
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed         0
eggmass.count <- eggmass.count[ !select, ] 


# Create a record for each species x each year
complete.species.year <- expand.grid(Study.Area.Name=unique(eggmass.count$Study.Area.Name),
                                     Transect.Label =unique(eggmass.count$Transect.Label),
                                     Species=unique(eggmass.count$Species),
                                     Year   =unique(eggmass.count$Year,  stringsAsFactors=FALSE))
eggmass.count <- merge(complete.species.year, eggmass.count, all.x=TRUE)

# replace all NA by 0
eggmass.count$n.eggmass[ is.na(eggmass.count$n.eggmass)] <- 0

# check that every species is given in each year
xtabs(~Species+Year, data=eggmass.count, exclude=NULL, na.action=na.pass)
##              Year
## Species       2013 2014 2015
##   A-AMGR         1    1    1
##   A-AMMA         1    1    1
##   A-ANBO         1    1    1
##   A-PSRE         1    1    1
##   A-RAAU         1    1    1
##   A-PRSE         1    1    1
##   ALL.species    1    1    1
# Section of code used for testing multiple transects#
#eggmass.count2 <- eggmass.count
#eggmass.count2$Transect.Label ='xx'
#eggmass.count2$n.eggmass <- rpois(1, eggmass.count$n.eggmass)
#eggmass.count2
#eggmass.count <- rbind(eggmass.count, eggmass.count2)

# Make a preliminary plot of total count by date
prelim.egg.plot <- ggplot(data=eggmass.count, aes(x=Year, y=n.eggmass, color=Transect.Label, linetype=Transect.Label))+
   ggtitle("Amphibian egg mass count data")+
   ylab("Amphibian egg masses")+
   geom_point(position=position_dodge(width=.5))+
   geom_line( position=position_dodge(width=.5))+
   facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
   scale_x_continuous(breaks=min(eggmass.count$Year,na.rm=TRUE):max(eggmass.count$Year,na.rm=TRUE))
prelim.egg.plot 

ggplot2::ggsave(plot=prelim.egg.plot, 
       file=paste(file.prefix,'-plot-prelim-egg.png',sep=""),
       h=6, w=6, units="in",dpi=300)


# This is a regression analysis with Year as the trend variable 
# In this case, there is only one transect, so it is not necessary to have the Transect.Label in the model
# We also don't need to model process (year specific effects) because there is only 1 transect

# Fit a linear regression for each species 
fits<- dlply(eggmass.count, "Species", function(eggmass.count){
   cat("\n\n\n *** Starting analysis for species ", as.character(eggmass.count$Species[1]), "\n")
   if(length(unique(eggmass.count$Transect.Label))>1){
      eggmass.count$Transect.LabelF <- factor(eggmass.count$Transect.Label)
      eggmass.count$YearF           <- factor(eggmass.count$Year)
      egg.fit <- lmerTest::lmer(log(n.eggmass) ~ Year + (1|Transect.LabelF) + (1|YearF) , data=eggmass.count)
      print(anova(egg.fit, ddf='kenward-roger'))
      print(summary(egg.fit))
      print(VarCorr(egg.fit))
   }
   if(length(unique(eggmass.count$Transect.Label))==1){
       # with only 1 transect, not necessary to put random effects effect in the model
       egg.fit <- glm(n.eggmass ~  Year , data=eggmass.count, family=quasipoisson)
       print(Anova(egg.fit, test="F", type=3))
       print(summary(egg.fit))
   }
   list(Species=eggmass.count$Species[1], fit=egg.fit)
})
## 
## 
## 
##  *** Starting analysis for species  A-AMGR 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##               SS Df      F Pr(>F)
## Year      1.5790  1 0.8759 0.5211
## Residuals 1.8028  1              
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##       1        2        3  
##  0.3296  -1.4631   0.6796  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1679.9760  1946.1905   0.863    0.547
## Year          -0.8341     0.9666  -0.863    0.547
## 
## (Dispersion parameter for quasipoisson family taken to be 1.802804)
## 
##     Null deviance: 4.2902  on 2  degrees of freedom
## Residual deviance: 2.7112  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## 
##  *** Starting analysis for species  A-AMMA 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##               SS Df          F    Pr(>F)    
## Year      4.3944  1 2.5127e+10 4.016e-06 ***
## Residuals 0.0000  1                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##         1          2          3  
##  2.11e-08  -1.87e-05  -2.11e-08  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 46621.82    2013.00   23.16   0.0275 *
## Year          -23.16       1.00  -23.16   0.0275 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 4.753932e-10)
## 
##     Null deviance: 4.3944e+00  on 2  degrees of freedom
## Residual deviance: 3.4978e-10  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 22
## 
## 
## 
## 
##  *** Starting analysis for species  A-ANBO 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##               SS Df          F    Pr(>F)    
## Year      2.1972  1 1.4968e+10 5.204e-06 ***
## Residuals 0.0000  1                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##          1           2           3  
## -2.110e-08  -1.713e-05   0.000e+00  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -45623.57    2015.00  -22.64   0.0281 *
## Year            22.64       1.00   22.64   0.0281 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 3.990352e-10)
## 
##     Null deviance: 2.1972e+00  on 2  degrees of freedom
## Residual deviance: 2.9359e-10  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 22
## 
## 
## 
## 
##  *** Starting analysis for species  A-PSRE 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##                SS Df      F Pr(>F)
## Year      0.51094  1 0.3558 0.6576
## Residuals 1.43614  1              
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##      1       2       3  
##  0.352  -1.353   0.555  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1052.1098  1823.2336   0.577    0.667
## Year          -0.5224     0.9054  -0.577    0.667
## 
## (Dispersion parameter for quasipoisson family taken to be 1.436144)
## 
##     Null deviance: 2.7726  on 2  degrees of freedom
## Residual deviance: 2.2616  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## 
##  *** Starting analysis for species  A-RAAU 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##            SS Df  F Pr(>F)
## Year      0.0  1  0      1
## Residuals 0.5  1          
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##       1        2        3  
## -0.3022   0.5372  -0.3022  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.877e-01  8.721e+02       0        1
## Year        1.050e-13  4.330e-01       0        1
## 
## (Dispersion parameter for quasipoisson family taken to be 0.5000005)
## 
##     Null deviance: 0.47113  on 2  degrees of freedom
## Residual deviance: 0.47113  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## 
##  *** Starting analysis for species  A-PRSE 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##           SS Df  F Pr(>F)
## Year       0  1  0      1
## Residuals 14  1          
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##      1       2       3  
## -2.160   2.459  -2.160  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  8.473e-01  3.488e+03       0        1
## Year        -4.488e-13  1.732e+00       0        1
## 
## (Dispersion parameter for quasipoisson family taken to be 14.00001)
## 
##     Null deviance: 15.381  on 2  degrees of freedom
## Residual deviance: 15.381  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## 
## 
##  *** Starting analysis for species  ALL.species 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.eggmass
## Error estimate based on Pearson residuals 
## 
##               SS Df      F Pr(>F)
## Year      1.1507  1 1.0721 0.4889
## Residuals 1.0733  1              
## 
## Call:
## glm(formula = n.eggmass ~ Year, family = quasipoisson, data = eggmass.count)
## 
## Deviance Residuals: 
##       1        2        3  
## -0.3711   0.8001  -0.5046  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 585.3617   569.3084   1.028    0.491
## Year         -0.2897     0.2827  -1.025    0.492
## 
## (Dispersion parameter for quasipoisson family taken to be 1.073323)
## 
##     Null deviance: 2.1832  on 2  degrees of freedom
## Residual deviance: 1.0325  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
# Get the model diagnostic plots
l_ply(fits, function(x){
    cat("\n\n\n *** Diagnostic plot for species ", as.character(x$Species), "\n")
   if(length(unique(x$fit$data$Transect.Label))>1){
      diag.plot <- sf.autoplot.lmer(x$fit)  # residual and other diagnostic plots
      plot(diag.plot)
   }
   if(length(unique(x$fit$data$Transect.Label))==1){
      diag.plot <- autoplot(x$fit)  # residual and other diagnostic plots
      show(diag.plot)
   }
   ggplot2:: ggsave(plot=diag.plot, 
          file=paste(file.prefix,"-egg-residual-plot-",x$Species,".png",sep=""),
          h=6, w=6, units="in", dpi=300)
})
## 
## 
## 
##  *** Diagnostic plot for species  A-AMGR

## 
## 
## 
##  *** Diagnostic plot for species  A-AMMA
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_segment).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).

## 
## 
## 
##  *** Diagnostic plot for species  A-ANBO
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).

## 
## 
## 
##  *** Diagnostic plot for species  A-PSRE

## 
## 
## 
##  *** Diagnostic plot for species  A-RAAU

## 
## 
## 
##  *** Diagnostic plot for species  A-PRSE

## 
## 
## 
##  *** Diagnostic plot for species  ALL.species

l_ply(fits, function(x){
   cat("\n\n\n *** Check for autocorrelation for species ", as.character(x$Species), "\n")
   eggmass.count <- x$fit$data
   # check for autocorrelation
   if(length(unique(eggmass.count$Transect.Label))>1){
      eggmass.count$resid <- log(eggmass.count$n.eggmass) - predict(x$fit, newdata=eggmass.count, re.form=~0)
   }
   #browser()
   if(length(unique(eggmass.count$Transect.Label))==1){
      eggmass.count$resid <- log(eggmass.count$n.eggmass+.1) - predict(x$fit, newdata=eggmass.count, type="link")
   }
   mean.resid <- plyr::ddply(eggmass.count, "Year", summarize, mean.resid=mean(resid))
   resid.fit <- lm( mean.resid ~ 1, data=mean.resid)
   dwres1 <- car::durbinWatsonTest(resid.fit)
   print(dwres1)
   dwres2 <- lmtest::dwtest(resid.fit)
   print(dwres2)
})
## 
## 
## 
##  *** Check for autocorrelation for species  A-AMGR 
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.6439684      2.931905   0.518
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 2.9319, p-value = 0.8818
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-AMMA 
##  lag Autocorrelation D-W Statistic p-value
##    1    -0.001098058      1.003294   0.122
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.0033, p-value = 0.01854
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-ANBO 
##  lag Autocorrelation D-W Statistic p-value
##    1   -0.0006940065      1.002082      NA
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.0021, p-value = 0.01282
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-PSRE 
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.6561264      2.968379   0.432
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 2.9684, p-value = 0.9197
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-RAAU 
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.6666667             3   0.166
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 3, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-PRSE 
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.6666667             3       0
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 3, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  ALL.species 
##  lag Autocorrelation D-W Statistic p-value
##    1       -0.642872      2.928616   0.444
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 2.9286, p-value = 0.8789
## alternative hypothesis: true autocorrelation is greater than 0
# extract a table of statistics for each study area x species combiations
egg.slopes <- ldply(fits, function(x){
  eggmass.count <- x$fit$data
  egg.fit       <- x$fit
  if(length(unique(eggmass.count$Transect.Label))==1){
     egg.slopes <- data.frame(
        Study.Area.Name =eggmass.count$Study.Area.Name[1],
        slope    = coef(egg.fit)["Year"],
        slope.se = sqrt(diag(vcov(egg.fit)))["Year"],
        p.value  = summary(egg.fit)$coefficients[row.names(summary(egg.fit)$coefficients)=="Year"  ,"Pr(>|t|)"],
        stringsAsFactors=FALSE
     )
   }
   if(length(unique(eggmass.count$Transect.Label))>1){
     egg.slopes <- data.frame(
          Study.Area.Name =eggmass.count$Study.Area.Name[1],
          slope           = fixef(egg.fit)["Year"],
          slope.se        = summary(egg.fit)$coefficients["Year","Pr(>|t|)"],
          p.value         = summary(egg.fit)$coefficients[row.names(summary(egg.fit)$coefficients)=="Year"  ,"Pr(>|t|)"], 
          stringsAsFactors=FALSE)
   }
   egg.slopes
})
egg.slopes
##       Species Study.Area.Name         slope  slope.se    p.value
## 1      A-AMGR Alice Lake Park -8.341152e-01 0.9665708 0.54674457
## 2      A-AMMA Alice Lake Park -2.316003e+01 0.9999993 0.02747079
## 3      A-ANBO Alice Lake Park  2.264197e+01 0.9999992 0.02809852
## 4      A-PSRE Alice Lake Park -5.224423e-01 0.9054296 0.66682897
## 5      A-RAAU Alice Lake Park  1.049719e-13 0.4330127 1.00000000
## 6      A-PRSE Alice Lake Park -4.487664e-13 1.7320508 1.00000000
## 7 ALL.species Alice Lake Park -2.896940e-01 0.2827022 0.49222412
# compute the fitted values from the model
egg.fitted <- ldply(fits, function(x){
   eggmass.count <- x$fit$data
   egg.fit       <- x$fit
   Transect.Label <- unique(as.character(eggmass.count$Transect.Label))
   newdata <- expand.grid(Year=seq(min(eggmass.count$Year, na.rm=TRUE),max(eggmass.count$Year, na.rm=TRUE), .1),
                          Transect.LabelF=Transect.Label)
   newdata$Transect.Label <-as.character(newdata$Transect.LabelF)
   if(length(unique(eggmass.count$Transect.Label))==1){
      newdata$pred.mean <- predict(egg.fit, newdata=newdata,type="response")
   }
   if(length(unique(eggmass.count$Transect.Label))>1){
      newdata$pred.mean <- exp(predict(egg.fit, newdata=newdata,type="response", re.form=~0))
   }
   # we now must average over all of the transect labels
   egg.fitted <- plyr::ddply(newdata, c("Year","Transect.Label"), plyr::summarize, pred.mean=mean(pred.mean))
   egg.fitted$Study.Area.Name <- eggmass.count$Study.Area.Name
   egg.fitted
})
egg.fitted
##         Species   Year     Transect.Label    pred.mean Study.Area.Name
## 1        A-AMGR 2013.0 Fawn Lake Transect 2.464816e+00 Alice Lake Park
## 2        A-AMGR 2013.1 Fawn Lake Transect 2.267563e+00 Alice Lake Park
## 3        A-AMGR 2013.2 Fawn Lake Transect 2.086096e+00 Alice Lake Park
## 4        A-AMGR 2013.3 Fawn Lake Transect 1.919151e+00 Alice Lake Park
## 5        A-AMGR 2013.4 Fawn Lake Transect 1.765566e+00 Alice Lake Park
## 6        A-AMGR 2013.5 Fawn Lake Transect 1.624272e+00 Alice Lake Park
## 7        A-AMGR 2013.6 Fawn Lake Transect 1.494285e+00 Alice Lake Park
## 8        A-AMGR 2013.7 Fawn Lake Transect 1.374702e+00 Alice Lake Park
## 9        A-AMGR 2013.8 Fawn Lake Transect 1.264688e+00 Alice Lake Park
## 10       A-AMGR 2013.9 Fawn Lake Transect 1.163478e+00 Alice Lake Park
## 11       A-AMGR 2014.0 Fawn Lake Transect 1.070368e+00 Alice Lake Park
## 12       A-AMGR 2014.1 Fawn Lake Transect 9.847087e-01 Alice Lake Park
## 13       A-AMGR 2014.2 Fawn Lake Transect 9.059049e-01 Alice Lake Park
## 14       A-AMGR 2014.3 Fawn Lake Transect 8.334076e-01 Alice Lake Park
## 15       A-AMGR 2014.4 Fawn Lake Transect 7.667120e-01 Alice Lake Park
## 16       A-AMGR 2014.5 Fawn Lake Transect 7.053540e-01 Alice Lake Park
## 17       A-AMGR 2014.6 Fawn Lake Transect 6.489062e-01 Alice Lake Park
## 18       A-AMGR 2014.7 Fawn Lake Transect 5.969759e-01 Alice Lake Park
## 19       A-AMGR 2014.8 Fawn Lake Transect 5.492014e-01 Alice Lake Park
## 20       A-AMGR 2014.9 Fawn Lake Transect 5.052502e-01 Alice Lake Park
## 21       A-AMGR 2015.0 Fawn Lake Transect 4.648162e-01 Alice Lake Park
## 22       A-AMMA 2013.0 Fawn Lake Transect 2.000000e+00 Alice Lake Park
## 23       A-AMMA 2013.1 Fawn Lake Transect 1.973344e-01 Alice Lake Park
## 24       A-AMMA 2013.2 Fawn Lake Transect 1.947044e-02 Alice Lake Park
## 25       A-AMMA 2013.3 Fawn Lake Transect 1.921094e-03 Alice Lake Park
## 26       A-AMMA 2013.4 Fawn Lake Transect 1.895490e-04 Alice Lake Park
## 27       A-AMMA 2013.5 Fawn Lake Transect 1.870227e-05 Alice Lake Park
## 28       A-AMMA 2013.6 Fawn Lake Transect 1.845301e-06 Alice Lake Park
## 29       A-AMMA 2013.7 Fawn Lake Transect 1.820707e-07 Alice Lake Park
## 30       A-AMMA 2013.8 Fawn Lake Transect 1.796441e-08 Alice Lake Park
## 31       A-AMMA 2013.9 Fawn Lake Transect 1.772498e-09 Alice Lake Park
## 32       A-AMMA 2014.0 Fawn Lake Transect 1.748875e-10 Alice Lake Park
## 33       A-AMMA 2014.1 Fawn Lake Transect 1.725566e-11 Alice Lake Park
## 34       A-AMMA 2014.2 Fawn Lake Transect 1.702568e-12 Alice Lake Park
## 35       A-AMMA 2014.3 Fawn Lake Transect 1.679876e-13 Alice Lake Park
## 36       A-AMMA 2014.4 Fawn Lake Transect 1.657487e-14 Alice Lake Park
## 37       A-AMMA 2014.5 Fawn Lake Transect 1.635397e-15 Alice Lake Park
## 38       A-AMMA 2014.6 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 39       A-AMMA 2014.7 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 40       A-AMMA 2014.8 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 41       A-AMMA 2014.9 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 42       A-AMMA 2015.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 43       A-ANBO 2013.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 44       A-ANBO 2013.1 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 45       A-ANBO 2013.2 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 46       A-ANBO 2013.3 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 47       A-ANBO 2013.4 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 48       A-ANBO 2013.5 Fawn Lake Transect 1.778588e-15 Alice Lake Park
## 49       A-ANBO 2013.6 Fawn Lake Transect 1.711606e-14 Alice Lake Park
## 50       A-ANBO 2013.7 Fawn Lake Transect 1.647146e-13 Alice Lake Park
## 51       A-ANBO 2013.8 Fawn Lake Transect 1.585113e-12 Alice Lake Park
## 52       A-ANBO 2013.9 Fawn Lake Transect 1.525417e-11 Alice Lake Park
## 53       A-ANBO 2014.0 Fawn Lake Transect 1.467969e-10 Alice Lake Park
## 54       A-ANBO 2014.1 Fawn Lake Transect 1.412685e-09 Alice Lake Park
## 55       A-ANBO 2014.2 Fawn Lake Transect 1.359482e-08 Alice Lake Park
## 56       A-ANBO 2014.3 Fawn Lake Transect 1.308283e-07 Alice Lake Park
## 57       A-ANBO 2014.4 Fawn Lake Transect 1.259013e-06 Alice Lake Park
## 58       A-ANBO 2014.5 Fawn Lake Transect 1.211598e-05 Alice Lake Park
## 59       A-ANBO 2014.6 Fawn Lake Transect 1.165968e-04 Alice Lake Park
## 60       A-ANBO 2014.7 Fawn Lake Transect 1.122057e-03 Alice Lake Park
## 61       A-ANBO 2014.8 Fawn Lake Transect 1.079800e-02 Alice Lake Park
## 62       A-ANBO 2014.9 Fawn Lake Transect 1.039134e-01 Alice Lake Park
## 63       A-ANBO 2015.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 64       A-PSRE 2013.0 Fawn Lake Transect 1.542573e+00 Alice Lake Park
## 65       A-PSRE 2013.1 Fawn Lake Transect 1.464051e+00 Alice Lake Park
## 66       A-PSRE 2013.2 Fawn Lake Transect 1.389527e+00 Alice Lake Park
## 67       A-PSRE 2013.3 Fawn Lake Transect 1.318796e+00 Alice Lake Park
## 68       A-PSRE 2013.4 Fawn Lake Transect 1.251665e+00 Alice Lake Park
## 69       A-PSRE 2013.5 Fawn Lake Transect 1.187952e+00 Alice Lake Park
## 70       A-PSRE 2013.6 Fawn Lake Transect 1.127481e+00 Alice Lake Park
## 71       A-PSRE 2013.7 Fawn Lake Transect 1.070089e+00 Alice Lake Park
## 72       A-PSRE 2013.8 Fawn Lake Transect 1.015619e+00 Alice Lake Park
## 73       A-PSRE 2013.9 Fawn Lake Transect 9.639206e-01 Alice Lake Park
## 74       A-PSRE 2014.0 Fawn Lake Transect 9.148542e-01 Alice Lake Park
## 75       A-PSRE 2014.1 Fawn Lake Transect 8.682854e-01 Alice Lake Park
## 76       A-PSRE 2014.2 Fawn Lake Transect 8.240871e-01 Alice Lake Park
## 77       A-PSRE 2014.3 Fawn Lake Transect 7.821387e-01 Alice Lake Park
## 78       A-PSRE 2014.4 Fawn Lake Transect 7.423255e-01 Alice Lake Park
## 79       A-PSRE 2014.5 Fawn Lake Transect 7.045389e-01 Alice Lake Park
## 80       A-PSRE 2014.6 Fawn Lake Transect 6.686758e-01 Alice Lake Park
## 81       A-PSRE 2014.7 Fawn Lake Transect 6.346382e-01 Alice Lake Park
## 82       A-PSRE 2014.8 Fawn Lake Transect 6.023333e-01 Alice Lake Park
## 83       A-PSRE 2014.9 Fawn Lake Transect 5.716727e-01 Alice Lake Park
## 84       A-PSRE 2015.0 Fawn Lake Transect 5.425729e-01 Alice Lake Park
## 85       A-RAAU 2013.0 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 86       A-RAAU 2013.1 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 87       A-RAAU 2013.2 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 88       A-RAAU 2013.3 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 89       A-RAAU 2013.4 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 90       A-RAAU 2013.5 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 91       A-RAAU 2013.6 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 92       A-RAAU 2013.7 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 93       A-RAAU 2013.8 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 94       A-RAAU 2013.9 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 95       A-RAAU 2014.0 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 96       A-RAAU 2014.1 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 97       A-RAAU 2014.2 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 98       A-RAAU 2014.3 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 99       A-RAAU 2014.4 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 100      A-RAAU 2014.5 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 101      A-RAAU 2014.6 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 102      A-RAAU 2014.7 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 103      A-RAAU 2014.8 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 104      A-RAAU 2014.9 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 105      A-RAAU 2015.0 Fawn Lake Transect 1.333333e+00 Alice Lake Park
## 106      A-PRSE 2013.0 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 107      A-PRSE 2013.1 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 108      A-PRSE 2013.2 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 109      A-PRSE 2013.3 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 110      A-PRSE 2013.4 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 111      A-PRSE 2013.5 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 112      A-PRSE 2013.6 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 113      A-PRSE 2013.7 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 114      A-PRSE 2013.8 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 115      A-PRSE 2013.9 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 116      A-PRSE 2014.0 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 117      A-PRSE 2014.1 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 118      A-PRSE 2014.2 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 119      A-PRSE 2014.3 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 120      A-PRSE 2014.4 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 121      A-PRSE 2014.5 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 122      A-PRSE 2014.6 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 123      A-PRSE 2014.7 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 124      A-PRSE 2014.8 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 125      A-PRSE 2014.9 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 126      A-PRSE 2015.0 Fawn Lake Transect 2.333333e+00 Alice Lake Park
## 127 ALL.species 2013.0 Fawn Lake Transect 9.095895e+00 Alice Lake Park
## 128 ALL.species 2013.1 Fawn Lake Transect 8.836173e+00 Alice Lake Park
## 129 ALL.species 2013.2 Fawn Lake Transect 8.583866e+00 Alice Lake Park
## 130 ALL.species 2013.3 Fawn Lake Transect 8.338764e+00 Alice Lake Park
## 131 ALL.species 2013.4 Fawn Lake Transect 8.100661e+00 Alice Lake Park
## 132 ALL.species 2013.5 Fawn Lake Transect 7.869356e+00 Alice Lake Park
## 133 ALL.species 2013.6 Fawn Lake Transect 7.644656e+00 Alice Lake Park
## 134 ALL.species 2013.7 Fawn Lake Transect 7.426372e+00 Alice Lake Park
## 135 ALL.species 2013.8 Fawn Lake Transect 7.214321e+00 Alice Lake Park
## 136 ALL.species 2013.9 Fawn Lake Transect 7.008324e+00 Alice Lake Park
## 137 ALL.species 2014.0 Fawn Lake Transect 6.808210e+00 Alice Lake Park
## 138 ALL.species 2014.1 Fawn Lake Transect 6.613810e+00 Alice Lake Park
## 139 ALL.species 2014.2 Fawn Lake Transect 6.424960e+00 Alice Lake Park
## 140 ALL.species 2014.3 Fawn Lake Transect 6.241503e+00 Alice Lake Park
## 141 ALL.species 2014.4 Fawn Lake Transect 6.063285e+00 Alice Lake Park
## 142 ALL.species 2014.5 Fawn Lake Transect 5.890155e+00 Alice Lake Park
## 143 ALL.species 2014.6 Fawn Lake Transect 5.721968e+00 Alice Lake Park
## 144 ALL.species 2014.7 Fawn Lake Transect 5.558584e+00 Alice Lake Park
## 145 ALL.species 2014.8 Fawn Lake Transect 5.399866e+00 Alice Lake Park
## 146 ALL.species 2014.9 Fawn Lake Transect 5.245679e+00 Alice Lake Park
## 147 ALL.species 2015.0 Fawn Lake Transect 5.095895e+00 Alice Lake Park
# Plot with trend line 
egg.plot.summary <- ggplot2::ggplot(data=eggmass.count,
                                    aes(x=Year, y=n.eggmass))+
   ggtitle("Amphibian eggmass count ")+
   ylab("Amphibian Eggmass Count")+
   geom_point(size=3, aes(color=Transect.Label))+
   facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
   geom_line(data=egg.fitted, aes(x=Year,y=pred.mean))+
   scale_x_continuous(breaks=min(eggmass.count$Year,na.rm=TRUE):max(eggmass.count$Year,na.rm=TRUE))+
   geom_text(data=egg.slopes, aes(x=min(eggmass.count$Year, na.rm=TRUE), y=max(eggmass.count$n.eggmass, na.rm=TRUE)), 
             label=paste("Slope (on log scale) : ",round(egg.slopes$slope,2), 
                         " ( SE "  ,round(egg.slopes$slope.se,2),")",
                         " \np :"    ,round(egg.slopes$p.value,3)),
                         hjust="left",vjust="top")+
   theme(legend.position="top")
egg.plot.summary

ggsave(plot=egg.plot.summary, 
       file=paste(file.prefix,'-egg-plot-summary.png',sep=""),
       h=6, w=6, units="in", dpi=300)
#  Analysis of the visual survey protocol

# Look at # adults 
# We need to assume that detection probabilities are consistent over time and space
# There is only one monitoring station at Alice Lake

dim(amph.df)
## [1] 38 26
amph.v.df <- amph.df[ amph.df$Survey.Type=='VI',]
dim(amph.v.df)
## [1] 35 26
# Count the total number of egg masses seen by species
adults.count.by.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
                       n.adults=sum(Life.Stage=='AD'))
adults.count.by.species
##    Study.Area.Name     Transect.Label Year     Species n.adults
## 1  Alice Lake Park Fawn Lake Transect 2013      A-AMGR        0
## 2  Alice Lake Park Fawn Lake Transect 2013      A-AMMA        0
## 3  Alice Lake Park Fawn Lake Transect 2013      A-ANBO        0
## 4  Alice Lake Park Fawn Lake Transect 2013      A-PSRE        0
## 5  Alice Lake Park Fawn Lake Transect 2013      A-RAAU        1
## 6  Alice Lake Park Fawn Lake Transect 2013   AMPHIBION        0
## 7  Alice Lake Park Fawn Lake Transect 2014      A-ANBO        1
## 8  Alice Lake Park Fawn Lake Transect 2014      A-PRSE        1
## 9  Alice Lake Park Fawn Lake Transect 2014      A-RAAU        0
## 10 Alice Lake Park Fawn Lake Transect 2015      A-AMGR        0
## 11 Alice Lake Park Fawn Lake Transect 2015      A-ANBO        3
## 12 Alice Lake Park Fawn Lake Transect 2015      A-PSRE        1
## 13 Alice Lake Park Fawn Lake Transect 2015      A-RAAU        0
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed        0
# Count the total number of egg masses seen for ALL species
adults.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year","Species"), summarize,
                       n.adults=sum(Life.Stage=='AD'))
adults.count.all.species <- plyr::ddply(amph.v.df, c("Study.Area.Name","Transect.Label","Year"), summarize,
                       n.adults=sum(Life.Stage=='AD'))
adults.count.all.species$Species <- 'ALL.species'

adults.count <- rbind(adults.count.by.species, adults.count.all.species)




# Need to impute 0 values for species not seen in a year


# check out records where adults count not identified
select <- adults.count$Species %in% c("AMPHIBION","Unidentifed") # Notice type in unidentified
adults.count[ select,]
##    Study.Area.Name     Transect.Label Year     Species n.adults
## 6  Alice Lake Park Fawn Lake Transect 2013   AMPHIBION        0
## 14 Alice Lake Park Fawn Lake Transect 2015 Unidentifed        0
adults.count <- adults.count[ !select, ] 


# Create a record for each species x each year
complete.species.year <- expand.grid(Study.Area.Name=unique(adults.count$Study.Area.Name),
                                     Transect.Label =unique(adults.count$Transect.Label),
                                     Species=unique(adults.count$Species),
                                     Year   =unique(adults.count$Year,  stringsAsFactors=FALSE))
adults.count <- merge(complete.species.year, adults.count, all.x=TRUE)

# replace all NA by 0
adults.count$n.adults[ is.na(adults.count$n.adults)] <- 0

# check that every species is given in each year
xtabs(~Species+Year, data=adults.count, exclude=NULL, na.action=na.pass)
##              Year
## Species       2013 2014 2015
##   A-AMGR         1    1    1
##   A-AMMA         1    1    1
##   A-ANBO         1    1    1
##   A-PSRE         1    1    1
##   A-RAAU         1    1    1
##   A-PRSE         1    1    1
##   ALL.species    1    1    1
# Section of code used for testing multiple transects#
#adults.count2 <- adults.count
#adults.count2$Transect.Label ='xx'
#adults.count2$n.adults <- rpois(1, 1+adults.count$n.adults)
#adults.count2
#adults.count <- rbind(adults.count, adults.count2)

# Make a preliminary plot of total count by date
# Make a preliminary plot of total count by date
prelim.adult.plot <- ggplot(data=adults.count, aes(x=Year, y=n.adults, color=Transect.Label, linetype=Transect.Label))+
   ggtitle("Amphibian adult count data")+
   ylab("Amphibian adults")+
   geom_point(position=position_dodge(width=.5))+
   geom_line( position=position_dodge(width=.5))+
   facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
   scale_x_continuous(breaks=min(adults.count$Year,na.rm=TRUE):max(adults.count$Year,na.rm=TRUE))
prelim.adult.plot 

ggplot2::ggsave(plot=prelim.adult.plot, 
       file=paste(file.prefix,'-plot-prelim-adult.png',sep=""),
       h=6, w=6, units="in",dpi=300)


# This is a regression analysis with Year as the trend variable 
# In this case, there is only one transect, so it is not necessary to have the Transect.Label in the model
# We also don't need to model process (year specific effects) because there is only 1 transect

# Fit a linear regression for each species 
fits<- dlply(adults.count, "Species", function(adults.count){
   cat("\n\n\n *** Starting analysis for species ", as.character(adults.count$Species[1]), "\n")
   if(length(unique(adults.count$Transect.Label))>1){
      adults.count$Transect.LabelF <- factor(adults.count$Transect.Label)
      adults.count$YearF           <- factor(adults.count$Year)
      adult.fit <- lmerTest::lmer(log(n.adults) ~ Year + (1|Transect.LabelF) + (1|YearF) , data=adults.count)
      print(anova(adult.fit, ddf='kenward-roger'))
      print(summary(adult.fit))
      print(VarCorr(adult.fit))
   }
   if(length(unique(adults.count$Transect.Label))==1){
       # with only 1 transect, not necessary to put random effects effect in the model
       adult.fit <- glm(n.adults ~  Year , data=adults.count, family=quasipoisson)
       print(Anova(adult.fit, test="F", type=3))
       print(summary(adult.fit))
   }
   list(Species=adults.count$Species[1], fit=adult.fit)
})
## 
## 
## 
##  *** Starting analysis for species  A-AMGR 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##                   SS Df  F Pr(>F)
## Year      0.0000e+00  1  0      1
## Residuals 2.2748e-10  1          
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
##          1           2           3  
## -1.231e-05  -1.231e-05  -1.231e-05  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -23.303   2466.636  -0.009    0.994
## Year           0.000      1.225   0.000    1.000
## 
## (Dispersion parameter for quasipoisson family taken to be 6.183461e-10)
## 
##     Null deviance: 0.0000e+00  on 2  degrees of freedom
## Residual deviance: 4.5495e-10  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 21
## 
## 
## 
## 
##  *** Starting analysis for species  A-AMMA 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##                   SS Df  F Pr(>F)
## Year      0.0000e+00  1  0      1
## Residuals 2.2748e-10  1          
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
##          1           2           3  
## -1.231e-05  -1.231e-05  -1.231e-05  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -23.303   2466.636  -0.009    0.994
## Year           0.000      1.225   0.000    1.000
## 
## (Dispersion parameter for quasipoisson family taken to be 6.183461e-10)
## 
##     Null deviance: 0.0000e+00  on 2  degrees of freedom
## Residual deviance: 4.5495e-10  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 21
## 
## 
## 
## 
##  *** Starting analysis for species  A-ANBO 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##               SS Df      F Pr(>F)
## Year      3.8586  1 13.083 0.1717
## Residuals 0.2949  1              
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
##        1         2         3  
## -0.55294   0.34401  -0.08681  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3048.0126  1064.8901  -2.862    0.214
## Year            1.5132     0.5285   2.863    0.214
## 
## (Dispersion parameter for quasipoisson family taken to be 0.294941)
## 
##     Null deviance: 4.29022  on 2  degrees of freedom
## Residual deviance: 0.43162  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## 
##  *** Starting analysis for species  A-PSRE 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##               SS Df          F    Pr(>F)    
## Year      2.1972  1 1.4968e+10 5.204e-06 ***
## Residuals 0.0000  1                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
##          1           2           3  
## -2.110e-08  -1.713e-05   0.000e+00  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -45623.57    2015.00  -22.64   0.0281 *
## Year            22.64       1.00   22.64   0.0281 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 3.990352e-10)
## 
##     Null deviance: 2.1972e+00  on 2  degrees of freedom
## Residual deviance: 2.9359e-10  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 22
## 
## 
## 
## 
##  *** Starting analysis for species  A-RAAU 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##               SS Df          F    Pr(>F)    
## Year      2.1972  1 1.4968e+10 5.204e-06 ***
## Residuals 0.0000  1                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
##          1           2           3  
##  0.000e+00  -1.713e-05  -2.110e-08  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) 45578.29    2013.00   22.64   0.0281 *
## Year          -22.64       1.00  -22.64   0.0281 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 3.990352e-10)
## 
##     Null deviance: 2.1972e+00  on 2  degrees of freedom
## Residual deviance: 2.9359e-10  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 22
## 
## 
## 
## 
##  *** Starting analysis for species  A-PRSE 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##           SS Df  F Pr(>F)
## Year       0  1  0      1
## Residuals  2  1          
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
##       1        2        3  
## -0.8165   0.9295  -0.8165  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.099e+00  3.488e+03       0        1
## Year         4.340e-13  1.732e+00       0        1
## 
## (Dispersion parameter for quasipoisson family taken to be 2.000103)
## 
##     Null deviance: 2.1972  on 2  degrees of freedom
## Residual deviance: 2.1972  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## 
##  *** Starting analysis for species  ALL.species 
## Analysis of Deviance Table (Type III tests)
## 
## Response: n.adults
## Error estimate based on Pearson residuals 
## 
##               SS Df          F    Pr(>F)    
## Year      2.0008  1 6.6857e+21 7.786e-12 ***
## Residuals 0.0000  1                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Call:
## glm(formula = n.adults ~ Year, family = quasipoisson, data = adults.count)
## 
## Deviance Residuals: 
## 1  2  3  
## 0  0  0  
## 
## Coefficients:
##               Estimate Std. Error    t value Pr(>|t|)    
## (Intercept) -1.395e+03  1.808e-08 -7.717e+10 8.25e-12 ***
## Year         6.931e-01  8.976e-12  7.722e+10 8.24e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 2.992646e-22)
## 
##     Null deviance:  2.0008e+00  on 2  degrees of freedom
## Residual deviance: -2.9926e-22  on 1  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 3
# Get the model diagnostic plots
l_ply(fits, function(x){
    cat("\n\n\n *** Diagnostic plot for species ", as.character(x$Species), "\n")
   if(length(unique(x$fit$data$Transect.Label))>1){
      diag.plot <- sf.autoplot.lmer(x$fit)  # residual and other diagnostic plots
      plot(diag.plot)
   }
   if(length(unique(x$fit$data$Transect.Label))==1){
      diag.plot <- autoplot(x$fit)  # residual and other diagnostic plots
      show(diag.plot)
   }
   ggplot2:: ggsave(plot=diag.plot, 
          file=paste(file.prefix,"-adult-residual-plot-",x$Species,".png",sep=""),
          h=6, w=6, units="in", dpi=300)
})
## 
## 
## 
##  *** Diagnostic plot for species  A-AMGR

## 
## 
## 
##  *** Diagnostic plot for species  A-AMMA

## 
## 
## 
##  *** Diagnostic plot for species  A-ANBO

## 
## 
## 
##  *** Diagnostic plot for species  A-PSRE
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).

## 
## 
## 
##  *** Diagnostic plot for species  A-RAAU
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_text).

## 
## 
## 
##  *** Diagnostic plot for species  A-PRSE

## 
## 
## 
##  *** Diagnostic plot for species  ALL.species

l_ply(fits, function(x){
   cat("\n\n\n *** Check for autocorrelation for species ", as.character(x$Species), "\n")
   adults.count <- x$fit$data
   # check for autocorrelation
   if(length(unique(adults.count$Transect.Label))>1){
      adults.count$resid <- log(adults.count$n.adults) - predict(x$fit, newdata=adults.count, re.form=~0)
   }
   #browser()
   if(length(unique(adults.count$Transect.Label))==1){
      adults.count$resid <- log(adults.count$n.adults+.1) - predict(x$fit, newdata=adults.count, type="link")
   }
   mean.resid <- plyr::ddply(adults.count, "Year", summarize, mean.resid=mean(resid))
   resid.fit <- lm( mean.resid ~ 1, data=mean.resid)
   dwres1 <- car::durbinWatsonTest(resid.fit)
   print(dwres1)
   dwres2 <- lmtest::dwtest(resid.fit)
   print(dwres2)
})
## 
## 
## 
##  *** Check for autocorrelation for species  A-AMGR
## Warning in summary.lm(model): essentially perfect fit: summary may be
## unreliable
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1666667           1.5      NA
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.5, p-value = 0.3328
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-AMMA
## Warning in summary.lm(model): essentially perfect fit: summary may be
## unreliable
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.1666667           1.5      NA
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.5, p-value = 0.3328
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-ANBO 
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.5254895      2.576468   0.426
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 2.5765, p-value = 0.6954
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-PSRE 
##  lag Autocorrelation D-W Statistic p-value
##    1   -0.0006940065      1.002082      NA
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.0021, p-value = 0.01282
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-RAAU 
##  lag Autocorrelation D-W Statistic p-value
##    1   -0.0006940065      1.002082    0.06
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.0021, p-value = 0.01282
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  A-PRSE 
##  lag Autocorrelation D-W Statistic p-value
##    1      -0.6666667             3      NA
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 3, p-value = 1
## alternative hypothesis: true autocorrelation is greater than 0
## 
## 
## 
## 
##  *** Check for autocorrelation for species  ALL.species 
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.02167572      1.065027   0.164
##  Alternative hypothesis: rho != 0
## 
##  Durbin-Watson test
## 
## data:  resid.fit
## DW = 1.065, p-value = 0.1138
## alternative hypothesis: true autocorrelation is greater than 0
# extract a table of statistics for each study area x species combiations
adult.slopes <- ldply(fits, function(x){
  adults.count <- x$fit$data
  adult.fit       <- x$fit
  if(length(unique(adults.count$Transect.Label))==1){
     adult.slopes <- data.frame(
        Study.Area.Name =adults.count$Study.Area.Name[1],
        slope    = coef(adult.fit)["Year"],
        slope.se = sqrt(diag(vcov(adult.fit)))["Year"],
        p.value  = summary(adult.fit)$coefficients[row.names(summary(adult.fit)$coefficients)=="Year"  ,"Pr(>|t|)"],
        stringsAsFactors=FALSE
     )
   }
   if(length(unique(adults.count$Transect.Label))>1){
     adult.slopes <- data.frame(
          Study.Area.Name =adults.count$Study.Area.Name[1],
          slope           = fixef(adult.fit)["Year"],
          slope.se        = summary(adult.fit)$coefficients["Year","Pr(>|t|)"],
          p.value         = summary(adult.fit)$coefficients[row.names(summary(adult.fit)$coefficients)=="Year"  ,"Pr(>|t|)"], 
          stringsAsFactors=FALSE)
   }
   adult.slopes
})
adult.slopes
##       Species Study.Area.Name         slope     slope.se      p.value
## 1      A-AMGR Alice Lake Park  0.000000e+00 1.224745e+00 1.000000e+00
## 2      A-AMMA Alice Lake Park  0.000000e+00 1.224745e+00 1.000000e+00
## 3      A-ANBO Alice Lake Park  1.513231e+00 5.285470e-01 2.139275e-01
## 4      A-PSRE Alice Lake Park  2.264197e+01 9.999992e-01 2.809852e-02
## 5      A-RAAU Alice Lake Park -2.264197e+01 9.999992e-01 2.809852e-02
## 6      A-PRSE Alice Lake Park  4.340178e-13 1.732051e+00 1.000000e+00
## 7 ALL.species Alice Lake Park  6.931472e-01 8.976133e-12 8.244113e-12
# compute the fitted values from the model
adult.fitted <- ldply(fits, function(x){
   adults.count <- x$fit$data
   adult.fit       <- x$fit
   Transect.Label <- unique(as.character(adults.count$Transect.Label))
   newdata <- expand.grid(Year=seq(min(adults.count$Year, na.rm=TRUE),max(adults.count$Year, na.rm=TRUE), .1),
                          Transect.LabelF=Transect.Label)
   newdata$Transect.Label <-as.character(newdata$Transect.LabelF)
   if(length(unique(adults.count$Transect.Label))==1){
      newdata$pred.mean <- predict(adult.fit, newdata=newdata,type="response")
   }
   if(length(unique(adults.count$Transect.Label))>1){
      newdata$pred.mean <- exp(predict(adult.fit, newdata=newdata,type="response", re.form=~0))
   }
   # we now must average over all of the transect labels
   adult.fitted <- plyr::ddply(newdata, c("Year","Transect.Label"), plyr::summarize, pred.mean=mean(pred.mean))
   adult.fitted$Study.Area.Name <- adults.count$Study.Area.Name
   adult.fitted
})
adult.fitted
##         Species   Year     Transect.Label    pred.mean Study.Area.Name
## 1        A-AMGR 2013.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 2        A-AMGR 2013.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 3        A-AMGR 2013.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 4        A-AMGR 2013.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 5        A-AMGR 2013.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 6        A-AMGR 2013.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 7        A-AMGR 2013.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 8        A-AMGR 2013.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 9        A-AMGR 2013.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 10       A-AMGR 2013.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 11       A-AMGR 2014.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 12       A-AMGR 2014.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 13       A-AMGR 2014.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 14       A-AMGR 2014.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 15       A-AMGR 2014.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 16       A-AMGR 2014.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 17       A-AMGR 2014.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 18       A-AMGR 2014.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 19       A-AMGR 2014.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 20       A-AMGR 2014.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 21       A-AMGR 2015.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 22       A-AMMA 2013.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 23       A-AMMA 2013.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 24       A-AMMA 2013.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 25       A-AMMA 2013.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 26       A-AMMA 2013.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 27       A-AMMA 2013.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 28       A-AMMA 2013.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 29       A-AMMA 2013.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 30       A-AMMA 2013.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 31       A-AMMA 2013.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 32       A-AMMA 2014.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 33       A-AMMA 2014.1 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 34       A-AMMA 2014.2 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 35       A-AMMA 2014.3 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 36       A-AMMA 2014.4 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 37       A-AMMA 2014.5 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 38       A-AMMA 2014.6 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 39       A-AMMA 2014.7 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 40       A-AMMA 2014.8 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 41       A-AMMA 2014.9 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 42       A-AMMA 2015.0 Fawn Lake Transect 7.582560e-11 Alice Lake Park
## 43       A-ANBO 2013.0 Fawn Lake Transect 1.528729e-01 Alice Lake Park
## 44       A-ANBO 2013.1 Fawn Lake Transect 1.778481e-01 Alice Lake Park
## 45       A-ANBO 2013.2 Fawn Lake Transect 2.069036e-01 Alice Lake Park
## 46       A-ANBO 2013.3 Fawn Lake Transect 2.407060e-01 Alice Lake Park
## 47       A-ANBO 2013.4 Fawn Lake Transect 2.800307e-01 Alice Lake Park
## 48       A-ANBO 2013.5 Fawn Lake Transect 3.257801e-01 Alice Lake Park
## 49       A-ANBO 2013.6 Fawn Lake Transect 3.790036e-01 Alice Lake Park
## 50       A-ANBO 2013.7 Fawn Lake Transect 4.409223e-01 Alice Lake Park
## 51       A-ANBO 2013.8 Fawn Lake Transect 5.129569e-01 Alice Lake Park
## 52       A-ANBO 2013.9 Fawn Lake Transect 5.967600e-01 Alice Lake Park
## 53       A-ANBO 2014.0 Fawn Lake Transect 6.942542e-01 Alice Lake Park
## 54       A-ANBO 2014.1 Fawn Lake Transect 8.076762e-01 Alice Lake Park
## 55       A-ANBO 2014.2 Fawn Lake Transect 9.396283e-01 Alice Lake Park
## 56       A-ANBO 2014.3 Fawn Lake Transect 1.093138e+00 Alice Lake Park
## 57       A-ANBO 2014.4 Fawn Lake Transect 1.271726e+00 Alice Lake Park
## 58       A-ANBO 2014.5 Fawn Lake Transect 1.479492e+00 Alice Lake Park
## 59       A-ANBO 2014.6 Fawn Lake Transect 1.721200e+00 Alice Lake Park
## 60       A-ANBO 2014.7 Fawn Lake Transect 2.002396e+00 Alice Lake Park
## 61       A-ANBO 2014.8 Fawn Lake Transect 2.329533e+00 Alice Lake Park
## 62       A-ANBO 2014.9 Fawn Lake Transect 2.710115e+00 Alice Lake Park
## 63       A-ANBO 2015.0 Fawn Lake Transect 3.152873e+00 Alice Lake Park
## 64       A-PSRE 2013.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 65       A-PSRE 2013.1 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 66       A-PSRE 2013.2 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 67       A-PSRE 2013.3 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 68       A-PSRE 2013.4 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 69       A-PSRE 2013.5 Fawn Lake Transect 1.778588e-15 Alice Lake Park
## 70       A-PSRE 2013.6 Fawn Lake Transect 1.711606e-14 Alice Lake Park
## 71       A-PSRE 2013.7 Fawn Lake Transect 1.647146e-13 Alice Lake Park
## 72       A-PSRE 2013.8 Fawn Lake Transect 1.585113e-12 Alice Lake Park
## 73       A-PSRE 2013.9 Fawn Lake Transect 1.525417e-11 Alice Lake Park
## 74       A-PSRE 2014.0 Fawn Lake Transect 1.467969e-10 Alice Lake Park
## 75       A-PSRE 2014.1 Fawn Lake Transect 1.412685e-09 Alice Lake Park
## 76       A-PSRE 2014.2 Fawn Lake Transect 1.359482e-08 Alice Lake Park
## 77       A-PSRE 2014.3 Fawn Lake Transect 1.308283e-07 Alice Lake Park
## 78       A-PSRE 2014.4 Fawn Lake Transect 1.259013e-06 Alice Lake Park
## 79       A-PSRE 2014.5 Fawn Lake Transect 1.211598e-05 Alice Lake Park
## 80       A-PSRE 2014.6 Fawn Lake Transect 1.165968e-04 Alice Lake Park
## 81       A-PSRE 2014.7 Fawn Lake Transect 1.122057e-03 Alice Lake Park
## 82       A-PSRE 2014.8 Fawn Lake Transect 1.079800e-02 Alice Lake Park
## 83       A-PSRE 2014.9 Fawn Lake Transect 1.039134e-01 Alice Lake Park
## 84       A-PSRE 2015.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 85       A-RAAU 2013.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 86       A-RAAU 2013.1 Fawn Lake Transect 1.039134e-01 Alice Lake Park
## 87       A-RAAU 2013.2 Fawn Lake Transect 1.079800e-02 Alice Lake Park
## 88       A-RAAU 2013.3 Fawn Lake Transect 1.122057e-03 Alice Lake Park
## 89       A-RAAU 2013.4 Fawn Lake Transect 1.165968e-04 Alice Lake Park
## 90       A-RAAU 2013.5 Fawn Lake Transect 1.211598e-05 Alice Lake Park
## 91       A-RAAU 2013.6 Fawn Lake Transect 1.259013e-06 Alice Lake Park
## 92       A-RAAU 2013.7 Fawn Lake Transect 1.308283e-07 Alice Lake Park
## 93       A-RAAU 2013.8 Fawn Lake Transect 1.359482e-08 Alice Lake Park
## 94       A-RAAU 2013.9 Fawn Lake Transect 1.412685e-09 Alice Lake Park
## 95       A-RAAU 2014.0 Fawn Lake Transect 1.467969e-10 Alice Lake Park
## 96       A-RAAU 2014.1 Fawn Lake Transect 1.525417e-11 Alice Lake Park
## 97       A-RAAU 2014.2 Fawn Lake Transect 1.585113e-12 Alice Lake Park
## 98       A-RAAU 2014.3 Fawn Lake Transect 1.647146e-13 Alice Lake Park
## 99       A-RAAU 2014.4 Fawn Lake Transect 1.711606e-14 Alice Lake Park
## 100      A-RAAU 2014.5 Fawn Lake Transect 1.778588e-15 Alice Lake Park
## 101      A-RAAU 2014.6 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 102      A-RAAU 2014.7 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 103      A-RAAU 2014.8 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 104      A-RAAU 2014.9 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 105      A-RAAU 2015.0 Fawn Lake Transect 2.220446e-16 Alice Lake Park
## 106      A-PRSE 2013.0 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 107      A-PRSE 2013.1 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 108      A-PRSE 2013.2 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 109      A-PRSE 2013.3 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 110      A-PRSE 2013.4 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 111      A-PRSE 2013.5 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 112      A-PRSE 2013.6 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 113      A-PRSE 2013.7 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 114      A-PRSE 2013.8 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 115      A-PRSE 2013.9 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 116      A-PRSE 2014.0 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 117      A-PRSE 2014.1 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 118      A-PRSE 2014.2 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 119      A-PRSE 2014.3 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 120      A-PRSE 2014.4 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 121      A-PRSE 2014.5 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 122      A-PRSE 2014.6 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 123      A-PRSE 2014.7 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 124      A-PRSE 2014.8 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 125      A-PRSE 2014.9 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 126      A-PRSE 2015.0 Fawn Lake Transect 3.333333e-01 Alice Lake Park
## 127 ALL.species 2013.0 Fawn Lake Transect 1.000000e+00 Alice Lake Park
## 128 ALL.species 2013.1 Fawn Lake Transect 1.071773e+00 Alice Lake Park
## 129 ALL.species 2013.2 Fawn Lake Transect 1.148698e+00 Alice Lake Park
## 130 ALL.species 2013.3 Fawn Lake Transect 1.231144e+00 Alice Lake Park
## 131 ALL.species 2013.4 Fawn Lake Transect 1.319508e+00 Alice Lake Park
## 132 ALL.species 2013.5 Fawn Lake Transect 1.414214e+00 Alice Lake Park
## 133 ALL.species 2013.6 Fawn Lake Transect 1.515717e+00 Alice Lake Park
## 134 ALL.species 2013.7 Fawn Lake Transect 1.624505e+00 Alice Lake Park
## 135 ALL.species 2013.8 Fawn Lake Transect 1.741101e+00 Alice Lake Park
## 136 ALL.species 2013.9 Fawn Lake Transect 1.866066e+00 Alice Lake Park
## 137 ALL.species 2014.0 Fawn Lake Transect 2.000000e+00 Alice Lake Park
## 138 ALL.species 2014.1 Fawn Lake Transect 2.143547e+00 Alice Lake Park
## 139 ALL.species 2014.2 Fawn Lake Transect 2.297397e+00 Alice Lake Park
## 140 ALL.species 2014.3 Fawn Lake Transect 2.462289e+00 Alice Lake Park
## 141 ALL.species 2014.4 Fawn Lake Transect 2.639016e+00 Alice Lake Park
## 142 ALL.species 2014.5 Fawn Lake Transect 2.828427e+00 Alice Lake Park
## 143 ALL.species 2014.6 Fawn Lake Transect 3.031433e+00 Alice Lake Park
## 144 ALL.species 2014.7 Fawn Lake Transect 3.249010e+00 Alice Lake Park
## 145 ALL.species 2014.8 Fawn Lake Transect 3.482202e+00 Alice Lake Park
## 146 ALL.species 2014.9 Fawn Lake Transect 3.732132e+00 Alice Lake Park
## 147 ALL.species 2015.0 Fawn Lake Transect 4.000000e+00 Alice Lake Park
# Plot with trend line 
adult.plot.summary <- ggplot2::ggplot(data=adults.count,
                                    aes(x=Year, y=n.adults))+
   ggtitle("Amphibian adults count ")+
   ylab("Amphibian adults Count")+
   geom_point(size=3, aes(color=Transect.Label))+
   facet_wrap(~interaction(Study.Area.Name,Species), ncol=2, scales="free_y")+
   geom_line(data=adult.fitted, aes(x=Year,y=pred.mean))+
   scale_x_continuous(breaks=min(adults.count$Year,na.rm=TRUE):max(adults.count$Year,na.rm=TRUE))+
   geom_text(data=adult.slopes, aes(x=min(adults.count$Year, na.rm=TRUE), y=max(adults.count$n.adults, na.rm=TRUE)), 
             label=paste("Slope (on log scale) : ",round(adult.slopes$slope,2), 
                         " ( SE "  ,round(adult.slopes$slope.se,2),")",
                         " \np :"    ,round(adult.slopes$p.value,3)),
                         hjust="left",vjust="top")+
   theme(legend.position="top")
adult.plot.summary

ggsave(plot=adult.plot.summary, 
       file=paste(file.prefix,'-adult-plot-summary.png',sep=""),
       h=6, w=6, units="in", dpi=300)